Electron-Electron Interactions in Isolated and Realistic Quantum Dots: 
A Density Functional Theory Study 
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We use Kohn-Sham spin-density-functional theory to study tiie statistics of ground-state spin and the spacing 
between conductance peaks in the Coulomb blockade regime for both 2D isolated and realistic quantum dots. 
We make a systematic investigation of the effects of electron-electron interaction strength and electron number 
on both the peak spacing and spin distributions. A direct comparison between the distributions from isolated and 
realistic dots shows that, despite the difference in the boundary conditions and confining potential, the statistical 
properties are qualitatively the same. Strong even/odd pairing in the peak spacing distribution is observed only 
in the weak e-e interaction regime and vanishes for moderate interactions. The probability of high spin ground 
states increases for stronger e-e interaction and seems to saturate around ~ 4. The saturated value is larger 
than previous theoretical predictions. Both spin and conductance peak spacing distributions show substantial 
variation as the electron number increases, not saturating until ~ 150. To interpret our numerical results, we 
analyze the spin distribution in the even TV case using a simple two-level model. 

PACS numbers: 73.23.Hk, 73.40.Gk, 73.63.Kv 
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I. INTRODUCTION 

The interplay of quantum mechanical interference and 
electron-electron interaction in semiconductor quantum dots 
(QDs) has attracted a lot of experimental and theoretical in- 
terest in recent years. Experimentally the ground state 
properties of a QD can be probed by Coulomb blockade (CB) 
measurements at near-zero temperature, in which the conduc- 
tance through the dot is blocked except at particular gate volt- 
ages, causing the conductance as a function of gate voltage 
to be a series of sharp peaks. In particular, the spacing be- 
tween neighboring peaks is related to the second difference 
of the ground state energy with respect to electron number 
N, A2E{N) = Eg,{N + l)+Eg,{N - l)-2Egs{N), which 
is often called the addition energy. The ground state spin of 
the dot can be inferred from the motion of the peak position 
in a magnetic field. ^'^ Except in small dots with high geo- 
metrical symmetry |i whose ground state properties show shell 
structure following Hund's rules, most experimental QDs are 
asymmetric due to either the irregular shape of the confining 
gate or impurity disorder, and the number of electrons is usu- 
ally several hundred . The properties of QDs show ap- 
parently random but reproducible fluctuations as some exter- 
nal parameters such as the chemical potential, the shape of the 
confining potential, external magnetic field and so on, vary. 
In this case, only statistical properties are of physical signifi- 
cance. 

While the statistical behavior of peak height fluctuations 
can be well described by a non-interacting random matrix 
theory (RMT) model,^ a similar model, the constant inter- 
action (CI) plus RMT, fails dramatically in describing the 
peak spacing distribution. The importance of e-e interac- 
tion was well recognized, but different theoretical approaches 
differ quantitatively or even qualitatively in many aspects 
of e-e interactions effects. The numerical studies for small 
dots^iiSii^Jiii^ii^iiii^ii^ yield a Gaussian-like peak spacing dis- 



tribution without any even/odd pairing effect, which agrees 
with experimental results on much larger dots. On the other 
hand, the statistical approaches for large dot£^°'^''^^ predict a 
much more pronounced even/odd effect at zero temperature 
that is absent in experimental results. 

Besides conductance peak spacings, another quantity that 
has attracted a lot of interest is the ground state spin in the 
CB regime. The occurrence of non-minimal ground state spin 
has been suggested as an explanation for the absence of a bi- 
modal distribution in conductance peak spacings,'^ and for 
the kinks in the parametric motion of the peaks. A recent 
experimental study on the change in CB peak position as a 
function of magnetic filed by Folk, et al.^ found that three 
out of five even-A^ states had S — 1. As in the case of peak 
spacing distributions, previous theoretical studies on ground 
state spin distributions in QDs gave very different results: 
Simulations on small tight-binding lattice models with ran- 
dom site energies predict a dominant fraction of high spin 
ground states, but results from KS-SDFT studies of small 
disordered parabolic quantum dots''' showed that the proba- 
bility of high-spin ground states never exceeds 50% even in 
the strong interaction regime. 

Despite intensive efforts to understand the role of e-e in- 
teraction in QD systems, there are still many questions that 
require further investigation. In this paper, we focus on the 
following issues: (1) The first one is concerned with the evolu- 
tion of statistical properties of QDs as the interaction strength 
increases. The e-e interaction strength is usually character- 
ized by the Wigner-Seitz radius r^, which is defined in 2D as 
rs = 1/ v^Tmao, where n is the electron density and gq is Bohr 
radius. Previous theoretical studies gave different answers to 
this question, indicating the necessity of further theoretical ex- 
ploration. In the meantime, this issue has attracted the experi- 
mental interest, and some relevant experiments are underway. 
(2) The second issue to be explored is the effect of confining 
potential. Previous theoretical studies were based on abstract 
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statistical assumptions, chaotic model external potentials, or 
Gaussian impurities. None of them have considered the realis- 
tic confining potential imposed in actual experiments. Though 
it is usually assumed that the statistical properties of QDs are 
of universal significance, the validity of that is still to be ver- 
ified. (3) The third issue that is worthy of further elaboration 
is how statistical properties of QDs evolve as the size of a QD 
system, mainly the electron number, increases, which is im- 
portant not only in solving the puzzling discrepancy between 
numerical small-dot calculations and semi-analytic large-dot 
statistical analysis, but also in understanding the fundamen- 
tal question in mesoscopic physics of how macroscopic prop- 
erties emerge from microscopic physics via the mesoscopic 
regime. 

In our previous study, ^^'^^ we reported a Kohn-Sham spin 
density functional theory (KS-SDFT) study of spin and con- 
ductance peak spacing distributions in 2D isolated quantum 
dots in a classically chaotic quartic oscillator potential with 
electron number N up to 200 at rg ^ 1.5. We found that for 
large electron number and asymmetric external potential, the 
peak spacing distribution is mainly Gaussian-like and there 
is no observable even/odd effects. We also found enhanced 
probability of high spin ground state implying stronger inter- 
action effects than expected from RPA approaches. In this pa- 
per, we extend the previous study in two directions: (1) Simi- 
lar calculations are done for weaker (r^ ~ 0.2) interactions to 
investigate how statistical properties of quantum dots evolve 
as the interaction strength varies. (2) We extended previous 
studies in isolated QDs, in which electrons are confined by 
an analytic chaotic potential without explicit consideration of 
the coupling with external gate, to QD systems in which the 
confining potential is calculated by solving the electrostatic 
Poisson equation with realistic boundary conditions. In the 
latter case, with more controllable parameters available, we 
are able to study statistical properties of QDs as a function of 
both e-e interaction strength (rg) and electron number (N) in 
a more systematic way. 

The outline of the paper is as follows: In the next section, 
the KS-SDFT method and the numerical techniques used in 
our study are concisely described. In Section III, we report the 
statistical properties of isolated quantum dots. In Section IV, 
the results for realistic quantum dots are presented. Section V 
concludes the paper with discussion about the significance of 
these results. 



II. METHOD 

A. Kohn-Sham spin density functional theory 

In KS-SDFT, the ground state energy of an interacting sys- 
tem with electron number N and total spin S in the local exter- 
nal potential V^xt (r) is written as a functional of spin densities 
n'^ with (T = a,(3 denoting spin-up and spin-down, respec- 
tively, 

E [n'',n^] = T, [n",n^] + / n(r)K>xt(r)dr 



n(r)n(r') 



(1) 



(Effective atomic units are used through the paper: for GaAs 
QDs with the dielectric constant e = 12.9 and effective elec- 
tron mass 771* — 0.06777ie, the values are 10.08 meV for en- 
ergy and 10.95 nm for length.) [n'',n'^] is the kinetic en- 
ergy of the KS non-interacting reference system which has 
the same ground state spin density as the interacting one, and 
Exc [n^jn^] is the exchange-correlation energy functional. 
The spin densities 77'^ satisfy the constraint J n'^{r)dr — N'^ 
with = {N + 2S) /2 and Nl^ ^ {N ~ 2S) /2. 

Assuming that the ground state of the non-interacting refer- 
ence system is non-degenerate, the non-interacting kinetic en- 
ergy is given by Tg [77", 77^] = J2i a {"^i |~^^^ | "^i)^ 
the ground state spin density is uniquely expressed as 



(2) 



Here ip" are the lowest single-particle orbitals which are ob- 
tained from 

{-iv2+V^,xt(r)+y^[77;r]+K-K,77^;r]}C(r) = e^CW, 

(3) 

where yff[77;r] and V^'^ [77" , 77^^ ; r] are the Hartree and 
exchange-correlation potentials, respectively. 



r ~ r 



(577'^(r) 



(4) 
(5) 



We have used the local spin-density approximation 
(LSDA)^^'^" for the exchange-correlation energy functional. 



i?XC ^",77''] 



?^(r)exc(?^(r),C(I■))rfI■, 



(6) 



where = (77" — 77^)/77 is the spin polarization, and exc is the 
exchange-correlation energy per electron in the homogeneous 
electron gas. We use Tanatar-Ceperley's parameterized form 
of Exc for the spin-compensated case (C, = 0) and fully spin 
polarized case {C, = 1).^' Exc for arbitrary spin polarization is 
obtained by the following interpolation,** 

excKC) =excKO)+/(C) [excKl) 



excKO)], (7) 



with 



/(C) = 



(i + C)'/^ + (i-C)'/^ 

23/2 „ 2 



(8) 



which is exact for the exchange energy density, but only 
approximate for the correlation energy density. Recently, 
Attaccalite, et al.^^ parameterized a new LSDA exchange- 
correlation functional form based on more accurate quantum 
Monte Carlo modeling data that includes spin polarization ex- 
plicitly. The new functional form was shown to be more ac- 
curate than Tanatar-Ceperley's in both zero and finite mag- 
netic field when compared to quantum Monte Carlo data.— 
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FIG. 1 : (Color online) Upper panel: Contours of Kxt in QOP system 
with parameters a = 10"*, 6 — 7r/4, A = 0.53, and 7 — 0.2. Lower 
panel: Electron density at TV = 200 and 5 = 0. It is obvious that the 
electron density in the isolated dot is far from uniform even within 
the classically allowed area; in the case shown here, the maximal 
density is 1 10% larger than that in the center. 



Numerical tests undertaken by us showed, however, that in 
the regime with which this paper is concerned, the difference 
between the two forms is minor In this paper, we will use 
Tanatar-Ceperley's form. Comparisons with quantum Monte- 
Carlo calculations for iV < 12 and interaction strength pa- 
rameter Ts < 8.0 have shown that LSDA works well for both 
the ground state spin and energy in 2D parabolic quantum 
dotsr^i^^^ as well as for many spin excitation energiesi^ 



B. Numerical techniques 

To investigate large electron number regime, we devel- 
oped efficient techniques for the implementation of KS-SDFT 
method which have been reported elsewhere. Essentially our 
algorithm includes the following components: (1) We use a 
particle-in-the-box basis for the representation of KS orbitals 
and for the action of the kinetic energy operator on the orbital 
wave functions; (2) The Kohn-Sham total energy is directly 
minimized by a conjugate gradient method that is based on 
Teter, Payne and Allan's method^^ but with important mod- 
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FIG. 2: (Color online) Illustration of data analysis procedure, (a) 
The addition energy A2E{N) as a function of TV (dotted), and its 
4"^ order polynomial fit that is related to increasing classical capac- 
itance as the number of electron in the dot increases (solid line), (b) 
The fluctuation of addition energy after removing the smooth trend, 
(c) Mean-level-spacing-scaled fluctuation of addition energy, whose 
statistical properties are to be studied. The data is from the QOP 
system with parameters a = 10"^, 6 = 7r/4, A = 0.53, and 7 = 0.2 



ifications, which include a more efficient line minimizaiton 
scheme and a delayed update of the effective potential;^^ (3) 
We use a Fourier convolution method to calculate Hartree po- 
tential efficiently;^^ (4) A simple one-way multi-grid tech- 
nique is used to accelerate the convergence.^'^ The high effi- 
ciency of our method makes the numerical study of statistical 
properties of large quantum dots with several hundred elec- 
trons become computationally feasiblei^iS^ 



III. SPIN AND PEAK SPACING DISTRIBUTIONS IN 
ISOLATED QUANTUM DOTS 

We first investigate statistical properties of quantum dots in 
a 2D isolated model system. Since experimental dots usually 
have irregular shapes that result in chaotic classical dynamics, 
we use the following quartic oscillator potential (QOP) as the 
external confining potential, 

K.t(r) =a[— + by^- 2Xx^y^ + j{x^y - xy')v\ . (9) 

Both the classical dynamics and the single-particle quantum 
mechanics at 7 = have been extensively studied:^" the sys- 
tem evolves continuously from integrable to fully chaotic as A 
changes from to 1 . The parameter 7 is introduced to break 
the four-fold symmetry. The prefactor a is used to adjust the 
flatness of the confining potential and hence the electron den- 
sity. 

Fig. ^shows the contours of a typical V^xt used in this study 
(upper panel), and the corresponding electron density at TV = 
200 (lower panel). A noticeable feature is that the electron 
density in isolated QD systems is not uniform: Electrons form 
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FIG. 3: (Color online) Distributions of dimensionless peak spacing [Eq. <10H for even (solid, red) and odd (dashed, blue) A'^ in isolated 
quantum dots at different interaction strengths. The inset of each figure shows the corresponding ground state spin distribution for A'^ even 
(solid, red) and odd (shaded, blue). A sliding window is used in estimating the probability density of peak spacing, yielding a smooth curve 
rather than a histogram; the width of the Gaussian window is taken asw— 5^^'',~^"""x where n = 200 is the number of data points. Note the 

° 2(l + Log2n) r 

striking evolution from clearly different even and odd distributions coupled with minimum ground state spin to very similar distributions with 
a large fraction of high spin ground states. 



peaks near the boundary region because of classical Coulomb 
repulsion between electrons. In real experimental systems, 
because of the coupling with the top gate, we expect to see 
an essentially uniform electron density distribution. Whether 
such a difference in electron density will result in significant 
deviation in statistical properties will be discussed in the next 
section. 

For a given Vcxt, the ground state energy Egg and spin 5gs 
as a function of N are determined by calculating several spin 
configurations for each N and selecting the one with minimal 
energy. The addition energy is then calculated as the second 
difference of EgsiN). As shown in Fig.Ha), A2E{N) de- 
creases slowly as N increases. This is mainly a classical effect 
- the increasing capacitance as the dot becomes large. Since 
we are mainly interested in the quantum mechanical effects, 
the classical smooth trend, denoted {A2E{N)), is fit by a 4- 
th order polynomial of N and then subtracted from A2E{N). 
To compare with experiments, the addition energy is further 
scaled by the mean level spacing. A, which is calculated by 
A — 2Trh'^ / m* Acff with Acff estimated as the classical al- 
lowed area; the resulting dimensionless spacing is denoted ^, 

. A.i.(A^) -^A2Em ^^^^ 

To obtain good statistics, we choose five sets of (A, 7) - (0.53, 
0.2), (0.565, 0.2), (0.6,0.1), (0.635, 0.15), and (0.67,0.1)- and 
for each set of parameters, we calculate a range of electron 
numbers N from 1 to 200. 

Fig. |3] plots peak spacing distributions in isolated QDs in 
different interaction regimes; ranges from about 0.2 to 1.3, 
and N — 120-200. In the large density or weak interaction 
regime, Ts ~0.2 [Fig.|3la)], the peak spacing distributions for 
even and odd N are sharply different: the distribution for even 
N is asymmetric and close to the Wigner surmise while that 
for odd is a narrow sharp symmetric peak, which agrees 
very well with the CI-RMT model. As the interaction strength 



increases, such an even/odd effect is gradually reduced: the 
odd N distribution becomes increasingly wider, and the even 
N distribution becomes more symmetric. At r., ~ 0.9 [Fig. 
|3jc)], the only difference between the two is a weak tail in 
the even N distribution at positive spacing. As is further 
increased to about 1.3 [Fig. |3jd)], which is the usual experi- 
mental regime, the even/odd effect is hardly perceptible, and 
the shape of the distribution is mainly Gaussian-like. 

The insets in Fig. I^show ground state spin distributions 
in the corresponding interaction regimes. At ~ 0.2, the 
minimal spin configuration dominates in both even and odd 
N cases, as one expects according to CI-RMT model. But 
even in this weak interaction regime, the portion of S" = 1 for 
even N is already significant, Ps=i =0.1. The probability of 
high-spin ground state increases quickly as the e-e interaction 
becomes stronger. At ^ 1.3, the high spin ground state 
is more probable than the minimal spin state. Such a result 
is remarkable considering that ^ 1.3 is only a moderate 
interaction strength. 



IV. SPIN AND PEAK SPACING DISTRIBUTIONS IN 2D 
REALISTIC QUANTUM DOTS 

A. Model 

The boundary condition used in isolated quantum dots is 
very different from that in the experimental systems. Fig. I^il- 
lustrates a typical experimental quantum dot structure, which 
is fabricated based on a GaAs-AlGaAs heterostructure. Nu- 
merical modeling of realistic quantum dotslislSiiiiil^ is still 
a great challenge in terms of computational efforts. The dif- 
ficulty mainly comes from the fact that although the volume 
where electrons are confined is small, the boundary conditions 
that couple the dot with the external confining gate are defined 



5 



f 

c 






f 

1 


4- 




s 


AlCJaAs 

- 




GaAs 




FIG. 4: (Color online) Schematic illustration of a realistic quan- 
tum dot studied in this paper. The dot is formed based on a GaAs- 
AlGaAs heterostructure, which consists of (from the bottom) an un- 
doped GaAs substrate, an undoped AlGaAs spacer layer, a n-doped 
AlGaAs layer, and finally a GaAs cap layer. The widths for the 
spacer, doped and cap layers are denoted as s, d and c, respectively. 
In this paper, s = 15 nm, d = 10 nm and c — 10 nm are taken. 




FIG. 5: Illustration of the shape of top confining gate used in this 
study. A negative voltage is imposed on the shaded region, depleting 
the electrons underneath so that the motion of electrons is confined 
to a small region. To obtain enough statistics, 20 different irregu- 
lar confining shapes are generated by changing xi, X2, yi and j/2 
randomly. 

on a much larger scale. Previous numerical studies of realis- 
tic quantum dots considered only small dots with less than a 
few tens of electrons, except in Ref. 43 where by reducing 
Schrodinger equations from 3D to multi-component 2D the 
author was able to study realistically a lateral quantum dot 
with N up to 100. 

We note that the most important factor in experimental con- 
ditions that is absent in the isolated quantum dot model is the 
coupling between the dot and the top gate. To take that into 
account, we introduce semi-infinite boundary conditions. "^''•'^^ 
Using R — (r. z) to denote the 3D spatial coordinate, R = 
(r, z) where z is the height in the growth direction and r the 
coordinates perpendicular to the growth direction, we take 



/)(r,z = 0) = Vg{r) + Vs.^ 
^1 = 



(11) 



where Vs is the Schottky barrier due to surface states, and 
Vg{r) is the gate voltage imposed on the top surface. In 2D 
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FIG. 6: (Color online) Upper panel: Contours of a sample 14, 
With the shape of the confining gate shown in Fig. |5| the confining 
potential is very irregular; its classical dynamics is believed to be 
chaotic. Lower panel: Electron density at N = 120. It is much more 
uniform that that in the isolated QOP shown in Fig. Q 

cases, the inclusion of semi-infinite realistic boundary condi- 
tions is quite straightforward. For Vcxt in Eq. Q, instead of 
an analytic form, it is calculated from 



V^^{v;z)^l- Jdr%{r') 



(|r-r'|2 + z2)3/2 



(12) 

where Vqw(z) is the confining potential in z-direction due to 
the quantum well or heterojunction structure from which the 
quantum dot is fabricated. In the case of 2D quantum dots, z 
is taken as the distance between the top surface and the mean 
electron height near the middle of the QD, zq. Because of the 
coupling with the gate, the Hartree potential has an additional 
image term^^ 



Vii{r) = / dr'n(r') 



1 



1 



r — r' 



-4z2)i/2J 



(13) 



This 2D quantum dot model will be denoted R2D to simplify 
the notation. 

Besides the more realistic boundary conditions included 
in the R2D model, its other advantage is that we can obtain 
statistics by changing the shape of the confining gate.'*'' In this 



FIG. 7: (Color online) Comparison of QOP and R2D models: dis- 
tributions of conductance peak spacing (solid for even TV, dashed for 
odd A'^) and ground state spin (insets) in (a) QOP at ~ 1.0; (b) 
R2D at Ts ~ 1.0; (c) QOP at Ts ~ 1.5; (d) R2D at ~ 1.5. For the 
QOP results = 80-160, while for the R2D model the additional 
gate averaging allows a smaller range A'^ = 110-130. Results for the 
two types of dots are in qualitative agreement. 



study, we use a confining gate with the shape illustrated in Fig. 
|5] By taking xi, X2, yi, and 2/2 randomly, we get a series of 
irregular confining potentials that are statistically independent 
as indicated by correlation analysis. Fig. |6lshows the contour 
of a typical confining potential and a sample electron density 
for ~ 200. We see that indeed the electron density in R2D 
is much more uniform than that in the isolated case. The avail- 
ability of a lot of statistically independent irregular potentials 
enables us to investigate how spin and peak spacing distribu- 
tions evolve as a function of N and in a more systematic 
way. 



B. Comparison of distributions from isolated and realistic 
quantum dots 

Fig. 0shows spin and peak spacing distributions from iso- 
lated and realistic QD systems at two values (r<; ^1.0 and 
1.5) within comparable electron number ranges. Qualitatively, 
they have similar features: At ~ 1.0, distributions from both 
QOP and R2D have a weak tail on the positive side for even 

compared to odd N, and at ^ 1.5, the tail vanishes. The 
spin distributions from two types of model systems are also 
the same within statistical uncertainty. Quantitatively there 
are indeed some observable differences in both spin and peak 
spacing distributions, but with the limited statistics available 
in this study, it is difficult to attribute physical significance to 
these subtle features. Our study therefore verifies the validity 
of using model systems to study statistical properties of quan- 
tum dots, and provides a direct demonstration of the universal 
feature of mesoscopic fluctuations. 



FIG. 8: (Color online) Probabilities of different spin states as a func- 
tion of rs for A*' = 11-40 in realistic quantum dots (circles: 5 = 0; 
squares: S — 1; diamonds: S = 2; up triangles: 5 = 1/2; down trian- 
gles: 5 = 3/2). The error bars shown in the plot are calculated from 
the standard bootstrap method. The ground state spin probabilities 
saturate at large Vs . 

C. Spin distributions as a function rs 

Fig. Oshows the spin distribution as a function of the inter- 
action strength in the electron number range N = 1 1-40. This 
study extends previous investigations for weak and moderate 
interaction strength in isolated dots to the stronger interaction 
regime. In agreement with previous studies, the fraction of 
high-spin ground states increases as the electron-electron in- 
teraction becomes stronger There are several features that are 
worth comment, (i) Crossing: For even N, the probability of 
5 = crosses with that of S* = 1 at ~ 3.0. (ii) Saturation: 
Probabilities of ground state spins for even N saturate at about 
rs ^ 3.5. The spin distribution for odd N is not saturated in 
the interaction strength range studied here, but the trend in- 
dicates that saturation is likely, (iii) Higher spin states: The 
probabilities of S = 2 and S — 5/2 remain small for all in- 
teraction strengths, though P{S — 2) becomes significant at 
rs > 3.0. 



D. Size Effects 

To investigate how statistical properties of quantum dots 
evolve as a function of electron number N, we calculated 
peak spacing and spin distributions in different electron num- 
ber ranges for N up to 200 at ~ 1.0 

Fig. |9] shows peak spacing distributions in different N 
ranges. For N — 11-40, the distributions for even and odd 
N are obviously different. The difference decreases as N in- 
creases, and becomes barely observable for N = 131-160 and 
N = 171-200. Our results show clearly that the conductance 
peak spacing distribution depends on electron number, and 
is significantly different between small {N < 50) and large 
(N > 100) dots. 

Fig. ^|shows the trend in ground state spin as N increases. 
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FIG. 9: (Color online) Distributions of dimensionless peak spacing for even (solid, red) and odd (dashed, blue) A'^ in realistic quantum dots 
in different A'^ ranges at rs ~ 1.0. There is a clear difference between even and odd distributions for small dots which, however, vanishes for 
large dots. 
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FIG. 10: (Color online) Probabilities of different spin states for dif- 
ferent A'^ ranges at Vs ~ 1.0 (circles: 5 = 0; squares: 5 = 1; up 
triangles: 5 = 1/2; down triangles: 5 = 3/2). The error bars shown 
in the plot are calculated from the standard bootstrap method. Note 
the substantial change from small to large dots. 

For this quantity, the dependence on N is even more signifi- 
cant: the probabiHty of high spin ground states increases as 
N increases. In particular, P{S = 1) increases from 0.27 at 
N = 11-40 to 0.45 at TV = 171-200. Another notable feature is 
that the spin distribution saturates at large N. 



E. Two-level model analysis 

KS-SDFT calculations reported above reveal some inter- 
esting features in conductance peak spacing and ground state 
spin distribution. To obtain a clearer picture of the underlying 
physics, and, in particular, to understand the main factors that 
determine the spin distribution in different regimes, we use a 
simple two-level model (TLM) to investigate the probability 
of singlet and triplet ground state spin for even N cases. The 
same model was used as an interpretive tool in Hirose and 
Wingreen's''^ work on spin distributions in small disordered 
parabolic quantum dots, and in doing that, they attribute satu- 
ration of the spin distribution to the effects of off-diagonal in- 



teraction matrix elements. Our analysis using single-particle 
wave functions in truly ballistic chaotic systems, however, 
leads to very different conclusions. 

The basic idea of the two-level model is to consider a two- 
electron system with only two single-particle levels. The in- 
teracting matrix elements (IME) 

V,,ki= J drdr'^*{r)rjir')v,,{r,r')Mr')Mr^)- (14) 

are calculated using a screened interaction kernel Wsc(r,r') 
to implicitly account for the effects of all other electrons. In 
particular, we use the screened interaction obtained from the 
2D random-phase approximation (RPA), 



Vsciq) = 



1 - Vciq)x{<i) 



(15) 



where Vc{q) = 27r/g is the Fourier transform of the bare 
Coulomb interaction, and x{q) is the susceptibility of the 2D 
electron gas. It is straightforward to include the screening due 
to the top gate by using Vc{q) = (27r/g)(l — c^^^"'^), but 
in the regime considered in this paper, the effects of the gate 
screening on the interacting matrix elements are negligible. 

In this TLM system, one can obtain singlet and triplet state 
energies either by exact diagonalization (ED) or in the frame- 
work of Hartree-Fock (HF) theory (ie., first-order perturbation 
theory). In TLM-HF, the difference of energy between singlet 
and triplet states has the following simple form. 



6E = Es=i - Es=o =Se- [Jn - ( J12 - K 



12) 



(16) 



Vijji are 



where Se is the single particle level spacing, Jij 
the direct Coulomb interaction terms, and K12 = V1212 is the 
exchange term. 

The main difference between our TLM analysis and Hi- 
rose and Wingreen's is in the way to choose single-particle 
energies and wave functions. In our analysis, we use single- 
particle eigen-energies and wavefunctions in the extended 
Thomas-Fermi effective potential^" that is obtained within the 
same realistic confining gate structure as in KS-SDFT calcu- 
lations. Wave functions obtained in this way are guaranteed 
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FIG. 1 1 : (Color online) Two-level model analysis of spin distribu- 
tion for TV = 1 1-40 as a function of interaction strength, (a) Probabil- 
ity of triplet ground state as a function of from TLM-ED (circles), 
TLM-HF (down triangles) and KS (up triangles), (b) Mean of diag- 
onal interaction matrix elements, Jn, J12, and K12, as well as the 
mean of the off-diagonal matrix element V1121, all in units of A. 

to be those of ballistic systems, and possible complications re- 
lated with the disorder strength (ie., localization) can therefore 
be avoided. We note that the two-level model with this way 
of calculating interacting matrix elements can be regarded as 
a simplified version of the model proposed by Ullmo, et al. in 
Ref. 20. 

Fig. Illf a) shows the probability of finding triplet ground 
states, P{S = 1), as a function of rg obtained from TLM-ED, 
TLM-HF, as well as KS-SDFT calculations. We see a striking 
agreement between TLM-ED and KS-SDFT results, which is 
indeed remarkable considering the extreme simplicity of our 
two-level model compared to KS-SDFT. Both TLM-HF and 
ED reproduce correctly the saturation of spin distribution at 
rs > 3.5. To explore possible mechanisms for the saturation, 
Fig.CUb) plots the mean values of Jn, J12, K12, and the off- 
diagonal term V1121, (the other off-diagonal term, V2221, is 
close to V1121, and is therefore not plotted here) as a function 
of rg. It is clearly shown that off-diagonal terms are an order 
of magnitude smaller than the Coulomb and exchange terms, 
indicating that the off-diagonal terms are negligible in terms 
of their contribution to the spin distribution. This is verified 
by the general agreement between TLM-HF and TLM-ED. 
On the other hand. Fig. II If b) shows that though ( J12) and 
{K12) do saturate at high r^, (Jn) does not. This indicates 
that considering mean values of Jn, J12, and K12 alone can 
not explain the saturation in the spin distribution. 

We also use the TLM to interpret the effect of finite size on 
the spin distributions, and the results are plotted in Fig. [T2I 
The simple two-level model also grasps most of the physics: 




80 120 160 

Electron Number 



200 



FIG. 12: (Color online) Two-level model (TLM) analysis of spin dis- 
tribution as a function of A'^. (a) Probability of triplet ground state as 
a function of A'^ from TLM-ED (circles), TLM-HF (down triangles) 
and KS (up triangles), (b) Means of diagonal interaction matrix ele- 
ments, Jii, J12, and K12, as well as that of the off-diagonal matrix 
element Vi 121, all in units of A. 

the probability of the triplet ground state calculated from both 
TLM-ED and HE increases as N increases, but saturates at 
large N, in agreement with the KS results. We note that the 
difference between different N ranges is mainly due to wave 
functions rather than eigenvalues (which are always random- 
matrix like for our system); for a ^/-electron system, the wave 
functions involved in the TLM are iJjn/2 and tpN/2+i- Our 
TLM results show clearly that the A^-dependence of the spin 
distributions comes mainly from the changing statistics of the 
wave functions. This is more directly demonstrated in Fig. 
I12r b). where mean values of Jn, J12, K12 and V^ii2i are plot- 
ted. We note that the A^-dependence of wave-function statis- 
tics is usually neglected in most statistical theories of quantum 
dots, where the large N limit is usually taken i^i^i^"i^^ Since 
many experimental data were obtained with N less than 200, 
our studies show clearly the necessity of considering finite N 
effects in interpreting those data. 



V. SUMMARY 

In this paper, we use Kohn-Sham spin density functional 
theory to study electron-electron interaction effects in quan- 
tum dots. In particular, we calculated conductance peak spac- 
ing and ground state spin distributions in an isolated quantum 
dot with classically chaotic quartic external potential and in a 
realistic quantum dot with explicit consideration of the cou- 
pling between the dot and the external confining gate. In the 
former case, as an extension of our previous study, we investi- 
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gated how the statistical properties evolve when the e-e inter- 
action strength increases from the weak (rg ~ 0.2) to moderate 
regimes (r^ 1.5) for N — 120-200. In the case of the realistic 
quantum dots, we first made a direct comparison between iso- 
lated and realistic quantum dots, and found that in spite of the 
difference in the boundary conditions, their statistical proper- 
ties are qualitatively similar. We further investigated spin and 
peak spacing distributions in strong interaction regimes with 
rg up to 4.0 and their dependence on the electron number N. 
The even/odd effect in peak spacing distributions vanishes as 
Tg increases and becomes weaker as N increases. 

We close with two comments on our results. First, our 
study, as well as other theoretical investigations, have made 
it quite clear that the lack of even/odd pairing in experimen- 
tal peak spacing distributions is due to the e-e interaction, and 
the CI-RMT model is therefore not valid for the description 
of mesoscopic fluctuations in Coulomb blockade peak spac- 
ings. In theoretical condensed matter physics, it is well ac- 
cepted that the e-e interaction at ^ 1.5 can be accurately 
described by many-body perturbation theory, most notably 
the RPA model. However, the results obtained from the RPA 
model combined with a RMT description of single particle 
properties disagree significantly from those obtained with nu- 
merical approaches. Understanding the sources for such dis- 



crepancies will deepen our general physical picture of quan- 
tum dots. 

Second, our results were obtained under the local spin den- 
sity approximation within KS-SDFT, whose validity, though 
well established in 3D bulk material systems through decades 
of experience, is not necessarily well justified for 2D large 
quantum dots with irregular confining potential and/or in the 
strong interaction regime. In the small dots for N up to about 
30, it is still computationally feasible to check the validity 
of LSDA by comparing with more accurate quantum Monte 
Carlo modeling, but not for the large-dot regime (N > 50). 
On the other hand, the statistical quantities investigated above 
are experimentally accessible. Comparing future experimen- 
tal data with LSDA predictions will provide a unique chance 
to check LSDA in real systems. 
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